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PCG SIGNAL ANALYZING 
USING WAVELET TRANSFORM 


S. GERGELY, ' M. V. PUSCA,”? M. N. ROMAN? 


Abstract. There are a very few PCG (Phonocardiography) signal characterization 
algorithms and also only a reduced number of devices capable for a complex pathology 
analysis. The analysis of the PCG signal is based on using of wavelet transform to 
achieve the full multi-resolution decomposition of the signal. The paper presents two such 
algorithms which are included ina portable device controlled by two microcontrollers. 
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1. Introduction 


Digital signal processing has a large number of analytical tools including the most 
important which is the Fourier analysis. Thus Fourier analysis is a technique that 
can transform a time domain signal into frequency domain. A number of 
applications in signal analysis require conservation or temporal information of a 
signal for which Fourier analysis encounters serious difficulties. To correct this 
deficiency Dennis Gabor introduced the concept of STFT (Short Time Fourier 
Transform) method which allows analysis of a shorter duration of the signal. This 
method solves in part the representation of the time-frequency signal, but there is 
a drawback by maintaining a constant time window throughout the frequency 
range for which the analysis is done. The appearance of wavelet transform allows 
the windows to be variable for the analyzed signal over the entire frequency 
spectrum. Therefore the first analyzing method displays a time-frequency 
scalogram of the PCG (Phonocardiography) signal. The second method enables 
the PCG signal analysis by converting two-dimensional convolution of a reference 
signal with the acquired signal [1]. By comparing the results of convolution it can 
be determined the correlation of each reference signal with the separately 
analyzed input signal. 


The basic functions used in the Fourier analysis are periodic functions (sine) with 
infinite duration. Thus it is possible to obtain information only in the frequency 
domain and of course the provided analyzed signal type must be periodically. 
Fourier transform for aperiodic signals will generate a continuous spectrum. 
Figure 1.1. presents a non-stationary signal which can be seen [6] that the inverse 
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Fourier transform permanently lose any information on the temporal position of 
the original signal. 


Time (s) 
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Frequency 


Time (s) 


Fig. 1.1. The Fourier transform of a aperiodic signal. 


By performing a transformation which uses a periodic basic function, well put out 
frequency spectral components but inevitably the time location information of the 
signal is lost. Also in the presence of [5] Gaussian noise in the analyzed signal, 
the isolation of time components becomes impossible because the present noise in 
the whole band cannot be eliminated by Fourier transform. In order to characterize 
a signal in time domain, it must be decomposed with aperiodic signals as a basis 
for well positioned time components. Therefore the wavelet transform is used in 
particular for analyzing signals that can be characterized as aperiodic type, 
containing noise, intermittent or having abrupt transitions. The strong ability of 
this type of transform [12] is that it can simultaneously examine both the time 
domain and frequency domain. This is done in a different way than the classical 
Fourier analysis. 
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Fig. 1.2. (a) Correlation of wavelet function with signal; 
(b) dilation-compression of wavelet function. 
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Either mentioned case can create a match of graphic form, compared with the 
original signal. If the wavelet function has a matching shape, then we get a great 
value for the correlation. On contrary if the wavelet function does not correlate 
with the portion of the signal then the correlation value of the transform decreases. 
From the mathematical point of view were studied a large number of types of 
wavelet functions but in scientific and engineering applications are used only up 
to four such functions. 


Of course it can be raised the question of why not to use the well-established 
Fourier methods for signal analysis. There are some major differences between 
Fourier analysis and wavelet transform performed by the method. Fourier analysis 
of the basic functions is localized in frequency but not in time. Local minor 
frequency changes are induced over the whole time domain of the Fourier 
transform. Wavelet functions are localized both in frequency and in time by 
shifting and scaling function without a projection over the entire transform. 
Another reason is that many functions classes (in our case signals) can be 
represented by this new method. Thus a signal having sharp discontinuities is 
using to approximate fewer wavelet functions than the sine-cosine functions of the 
typical Fourier transform. 


2. The wavelet transform 


The fundamental characteristic of wavelet functions is their ability to deal with 
noisy signals even covered by a lot of information, because in the end after the 
transform, the entire signal is compressed only to the wavelet coefficients only. 
Another important aspect is that in terms of number of necessary calculations for 
the Fourier transform there are required O(n-log,(n)) operations, instead of O(n) 


operations required by the wavelet transform. To complete the comparison 
between Fourier and wavelet transform should be noted that once the wavelet 
function y(r) is defined it can take the translation and _ the 


dilation (=?) e0)- R* “Rh. By defining the wavelet basis, a and b will take 
a 


the values of a=27/ and b=k-2/ where k andj eR. 


To fit a function in the class of wavelet type, they must meet the following 
conditions [15]: 


a) The function must have finite energy: 
E= [W(t}'dr<x (2.1) 


If the function y(t) is complex then the amplitude is computed using both real and 
imaginary parts. 
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b) if AG ) is the Fourier transform of function y(t) then: 
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And also the condition: 
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Where the above condition is the admisibility function which implies that the 


AN 
function has no zero frequency components y(0)=0 or the average value of y/(r) 
must be equal zero. It turned to be obvious that each candidate to a wavelet based 
class will have an admissibility constant Cg. 


For a function belonging to the space of finite energy y </’ (R) we get 


the transform: 


= (2.4) 


In order to ensure that at each scale factor "a" the wavelet function will have the 
same energy, energy conservation w(a) introduces a coefficient with value 1/ (PA 


Analog to the Fourier analysis, wavelet transform has the ability to reconstruct the 
original signal by computing the inverse wavelet transform of the analyzed signal. 
The mathematical form of the inverse wavelet transform is described by equation 


the signal: 


1 rr dadb 
a(t)= at [[7@)-v.0 3 
a 

% =00 (2.5) 
where the C, coefficient defines the admissibility condition. If we choose a narrower 
area of integration over the reverse wavelet transform it is obtained a filtering effect 
of the original signal. The above equation is called also a convolution therefore the 
final form of the wavelet transform is given by the following equation. This is also 
the well-known equation used in all practical applications: 


Vaplt)= su) a 
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3. Implementation of the filtering algorithm 


Stephane Mallat introduces the concept of filtering a signal through the discrete 
wavelet transformed in a multi-resolution decomposition. Each iteration of the 
process is called octave that can be assimilated to a set of two FIR filter type. By 
direct application of wavelet transform on a discrete signal, we obtain a double 
number of samples for the half of the transform vector. Because the number of 
points after each transform iteration is.equal to the number of samples of the 
original signal, we had to proceed to decimation by two.of the resulted vector. 


The procedure of signal decomposition is the reverse process of synthesis by 
wavelet transform. This is done by inserting the inverse low pass filters 
respectively high pass filters. A low-pass filter generates a mean signal (of two 
adjacent samples), while a high pass filter generates the detail elements by making 
the difference between two adjacent samples [1]. If case of a low-pass filter with 
coefficients {I/2,1/2}, the result is (x[n]+x{n—-1))/2, this obviously is the average of 
two samples. A high-pass filter having the {I/2,-1/2} yields the result of 
(x{n]— x{n—1])/2. The multi-resolution analysis introduces the average of a signal in 
a different set of filters as shown in figure 3.1. 


Decompozition Reconstruction 


= 2. 


(b) Wavelet synthesis and reconstruction algorithm 


Fig. 3.1. (a) Decomposition algorithm (b) multi-resolution synthesis algorithm 


The above described algorithm is part of a filtering process done by decimating 
the signal by two in which a sub-sampling is accomplished by removing 
sequentially one sample of the signal. The detail coefficients are obtained in a 
similar manner by using the pre-computed high pass filtering coefficients. The 
scaling sequence has a low pass filter behavior and the wavelet sequence is a low 
pass filter both where both are normalized for a coefficient sum equal with 2. In 
some applications the filtering sequences are normalized to have the sum equal 


with J2 , to remain orthogonally to each other. 
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4. Choosing the optimal wavelet function for the analysis of PGC signals 


Every application which requires the using of a wavelet transform requires the 
adaptation or finding of a certain type of wavelet function. The first criterion is 
defined by the correlation factor between the wavelet function and the analyzed 
signal. Unfortunately the shape of the correlated wavelet function is capable to 
indicate only the family of wavelet functions but gives no indication regarding the 
number of the wavelet coefficients...To..compute the necessary coefficient 
numbers, we used the wavelet toolbox of the Matlab program. A key criterion in 
determining the type of wavelet function consists in evaluation the overall 
reconstruction error, for the given type of signals that will appear in the 
application. So the evaluation must be made for signal reconstruction at the 
extreme in terms of spectral content. We imposed this condition for the 
decomposition algorithm to use a minimum number of coefficients, to optimize 
the convolution calculation time. Figure 4.1 is presenting the minimal set 
necessary for the calculation of the reconstruction error. 


load Normal; %Load file 
s=Normal(1:32768); 
[C,L]=wavedec(s,3,'db2') 
A3=wrcoef('a',C,L,'db2',3 
D1=wrcoef('d',C,L,'db2',1 
D2=wrcoef('d',C,L,'db2',2 
D3=wrcoef('d',C,L,'db2',3 
sig=A3+D3+D2+D1; %Reconstruction from signal parts 
err=max(abs(s-sig)); %Error computation 


%Computing the decomposition of level 3 


) 
) 
) 
) 


Fig. 4.1. Matlab program for wavelet reconstruction error computing. 


After computing the wavelet function we tested the validity of the selection by taking 
on all detail components of reconstructed signal and computing the reconstruction 
error. The goal was to find a maximum error therefore the maximum sensitivity of the 
wavelet function to the given input signal. The tables are presenting the overall 
computed errors which must be the lowest possible and the component extraction 
error which must indicate the highest sensitivity of the wavelet functions. 


These two conditions are met by the Daubeschies db2 function. It seemed that this 
function did answer very well to the PCG signal which is a strongly polynomial 
type due the massive presence of arbitrary heart murmur signals. Due to the 
Matlab convention the db2 wavelet function is mathematically defined by four 
coefficients which are computed by using the orthogonally conditions along with 
a specially imposed condition which gives the smoothness of the Daubeschies 
functions. The signal is reconstructed by using the equation: 

S = A(level ) + D(level ) + D(level _ 1) + D( level — 2) 


reconstrut 


(4.1) 


96 S. Gergely, M. V. Pusca, M. N. Roman 


Table 4.1. Reconstruction error computation for a normal PCG signal (without pathology) 


OVERALL 
DAUBESCHIES AVERAGE 
WAVELETS S=A3+D3+D2 | S=A3+D3+D1 | S=A3+D1+D2 ERROR ERROR 
S=A3+D3+D2+D1 


DB2 | DB2 (MATLAB) | 0. | 0.095 | 0.0168 0.0688 0.0317 1.46E-13 


Table 4.2. Reconstruction error computation for a PCG signal with Aortic regurgitation 


OVERALL 
ates. | S=A3+D3+D2 | S=A3+D3+D1 | S=A3+D1+D2 | AYERAGE ERROR 
S=A3+D3+D2+D1 


DB2 (MATLAB) 0.0678 0.1796 0.5198 0.255 1.168E-12 
DB3 (MATLAB) 0.035 0.151 0.5226 0.2362 1.116E-11 


| DB4 (MATLAB) | | DB4 (MATLAB) | | 0.0408 | 0408 | 0.1015 | 1015 0. os, | eae 171 De | -2.068E-12 | 12 


hate) 0.0596 0.0761 0.3775 0.171 3.12E-12 


| DB6 (MATLAB) | | DB6 (MATLAB) | 0. | 0.0646 0. | 0.0869 y= 0389S 3311 1. | 1.793E-12—_ | 12 


Since the reconstruction is made more difficult for ae with a higher 
spectral content so the finally chosen wavelet function is the db2. The reconstruction 
of the initially decomposed signal by using the wavelet transform is done through the 
reconstruction filters. The number of coefficients which are defining the wavelet 
function has a major influence over the construction of the reconstruction filter. 
According to Strang and Nguyen the reconstruction filter presented in figure 2.10, is a 
quadrature filter used very frequently in signal processing. Decimating the signal with 
a factor of two induces a series of aliased samples in the original signal which had to 
be eliminated through a anti-alias filter before the sampling process. The 
reconstruction signal [7] in terms of levels of decomposition is shown in figure 2.15. 


S=A1+D1 Level 1 


S=A2+D2+D1 Level 2 


S=A3+D3+D2+D1 Level 3 


Fig. 2.15. Reconstruction of signal in terms of decomposition levels 
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Through the PCG signal analysis, ultimately the reconstruction algorithm is used 
only in the process of choosing the right wavelet function for the multi-resolution 
decomposition. 


5. Algorithms for the characterization of the PCG signals 


Due to the extremely complex structure of the PCG signal, in analogy with an 
information carrier signal, it is characterized by a low frequency amplitude 
modulation, followed by a frequency modulation of a "carrier" generated by the 
cardiac activity at different times in relation to the movement of valves, and 
finally the overall signal is phase modulated between all components. The 
properties of the PCG signal in the frequency domain are the best starting point in 
the characterization of these signals. The most important parameter is the 
frequency domain of the signal which is situated in the range of 62-800 Hz. 
Figures 5.1 and 5.2 are showing the frequency spectrum domain for two cases. 


Oo 100 200 300 400 500 600 700 800 900 
Frequency 


Fig. 5.1. Normal heart PCG signal frequency spectrum. 


0 100 20 30 40 500 600 700 880 SO 1000 
Frequency 


Fig. 5.2. Aortic Regurgitation PCG signal frequency spectrum. 
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Figure 5.1 and 5.2 are showing that the frequency spectrum for the PCG signals 
are situated in the low and in the middle of the 62-800 Hz domain, regardless the 
type of the pathology. So the presentations of analysis algorithms of the PCG 
signals are presented through the following steps: 

e determination of the spectral content of the wavelet packets 

e precise detection of the heart beat rate (HBR) for the signals containing 

certain pathologies 

e PCG signal rescaling through multirate sampling 

e PCG signal envelope detection 

e Correlation of the two dimensional converted PCG signals 


Decomposition level Algorithm index 


Energy value of samples 
S*=) ms | | ks* =! 
0,0 Address 


Fig. 5.3. Computing algorithm of the time-frequency scalogram display addreses. 


Displaying of the time-frequency scalogram requires the construction of a special 
algorithm for the computing of addresses for each pixels of the color LCD 
display. The origin of the first displayed value is situated in the upper left corner 
of the display panel. The mentioned starting point is not appropriate for the 
display of the scalogram due the mathematical configuration of the time- 
frequency representation. The scalogram changes the standard displaying mode to 
a custom designed algorithm which modifies both horizontal (time) and vertical 
(frequency) coordinates. The main issue is the correct representation of the 
frequency coordinates, due to the interlaced addressing of the pixels. This 
addressing mode is imposed by the changing of the usual order of the computed 
frequency domains. The changing of the domain order follows a Gray code 
pattern; therefore each Gray code of the address value is first converted to its 
decimal value. The new value represents the address of the real pixel rows of the 
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time-frequency scalogram. The LCD display has 128 horizontal pixels and 160 
vertical pixels. Therefore finally the y coordinate yields by computing 
INDEX adressé = 159 - valeray: 


The value written at a specific address is given by the energy of the aquired 
samples, where the energy is given by the squared value of each filtered sample. 
The energy values are scaled to different colors which are represented in a color 
matrix coded in a RGB565 configuration. 


The scalogram algorithm is presented in figure 5.3. 


As a result of the full decomposition, the displayable resolution is given by 
equation 5.1: 
ee AM and _ 800-62 
nr[INDEX,,.| 128 


=5.76Hz 
(5.1) 


By adding new levels to the decomposition obviously the frequency resolution 
does increase by adding more frequency bands but also the time resolution 
decreases due the decreasing of the samples number. This is a practical prove of 
the Heisenberg principle of uncertainty which states that a time-frequency 
representation has a fundamental limitation because the product of time resolution 
At and frequency resolution Af is given by the equation 5.2: 


At-A@> i 
2 (5.2) 
Where: 
[ewofa 
OY TvoFar 
t 
Jv (5.3) 
fo ¥(of do 

Jy (5.4) 


and where ‘(@) is the Fourier transform of basis function (7). 


Therefore Heisenberg's uncertainty principle states that the frequency resolution 
can be achieved at the expense of resolution in time and vice versa. This is 
because the resolutions cannot be arbitrarily small because their product has a 
finite value. In figure 5.4 are presented some scalograms achieved by simulating 
the decomposition algorithm presented above. Simulations have shown the need 
for a suitable scaling corresponding to the color values assigned to the displayed 
energy. 
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Normal heart rate Atrial fibrilation 
Pulmonary stenosis Aortic stenosis Friction of the pericard 


All images are displayed on 128x128 pixels 


Fig. 5.4. Time-frequency scalograms obtained using algorithm presented in figure 5.3. 


The final form of the PCG signal correlation does include some intermediary 
algorithms presented in figure 5.5. 


: |= | Time translation of 
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parameter 


Correlation 
of PCG 
signal 


Fig. 5.5. PCG signal correlation algorithm. 
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Due to the large overall mathematical complexity, the algorithms are only 
enumerated in this paper. The above presented scalogram and the next algorithms 
are by now implemented in a portable medical device capable to characterize a 
number of pathologies which are present in the PCG signals. 


Simulations showed that the classical correlation method between a one- 
dimensional signal vectors with another reference signal, did not give the best 
results. The problems were caused by.the presence of relative symmetry of the 
shape in the reference signal. For example, even if the reference signal contains 
two asymmetric peaks with values which are very different, the same set of values 
can be found in any other signal due to the summing operation of the correlation 
process. 


Due to the complex nature of the PCG signal characterized by the presence of 
several peaks situated at relatively the same distance but in different pathology, it 
was imposed the introduction of an additional information which to indicate the 
relative position of signal components. Dividing the signal into successive epochs 
by applying threshold values and classical comparing applied to ECG signals can 
not be used. Tests have shown that signal synchronization using threshold method 
can be successful only when there's no pathology presence, and only to determine 
heart rate. Accurate and nondestructive detection by filtering the envelope of the 
PCG signal, pave the way for much closer analysis to a visual perception so that 
formal representation can be similar to an ECG signal. This requires the obtaining 
of separate images for each type of pathology stored in the SD card. 


HL 


Normal heart reference (N Mitral regurgitation reference (MR) S4 reference 
Aortic stenosis reference (AS) Mitral stenosis reference (MS) Pulmonary stenosis reference (PS) 


Fig. 5.6. PCG signal pathologies references. 
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Defining certain pathologies is done by manually extracting the desired signal 
sequence. The image extraction has to be done by a medically trained person who 
knows the components that make up a given pathology. The images were 
extracted from the demodulated signals, for which it was wrote an own program. 
Each reference image is extracted for a specific cardiac pathology. This can build 
a database with various pathologies to be compressed and stored. 


The necessary memory to store a pathology image is approximately 2 Ko. For a 
correct reference image generation, the computer requires as input parameter only 
the value of heart rate taken at time of the signal acquisition. The program extracts 
from the PCG signal a number of samples by moving two cursors on the signal 
and after that it allows to save as a file in a reference data or graphical bmp 
format. Some of the images corresponding to the pathologies references are 
shown in figure 5.6. 


During the researches we applied a new method to analyze the PCG signal by 
introducing an additional coordinate to the signal vector. The idea is that we 
switched to two-dimensional analysis similar to pattern recognition located in an 
image. We converted from one-dimensional vector v [x] to the vector w[y][x] by 
marking the position (indexes) sample with value 1, if the signal from that 
amplitude value is in the rescaled amplitude domain. Samples that fall outside the 
scope of amplitude are marked with 0. Thus the points defined in terms of 
measured values XOY are introducing a second dimension of the signal and can 
thus be interpreted as an image. Reducing the number of calculations required to 
scale the amplitude vector x [t] for a number of dots (pixels). By extending the 
method, if the marked pixel value changes from | to higher value, the vector 
wlLy][x], receives the z coordinate, i.e. the intensity that we refer to a grayscale or 
color image. 


Amplitude Signal reference 


Sample number [x] 


v[10]={1, 2, 4, 6, 8, 9, 7, 3, 1, 0} 


Fig. 5.7. Two-dimensional conversion of the PCG signal. 
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This type of resulted image can be processed by using advanced image processing 
methods. On a closer analysis, there is an important feature which allows to limit 
one of the coordinates of the vector w[y][x]. It can be observed that both signals 
involved in the algorithm are scaled to the same number of values on the Y 
coordinate, which provides a significantly higher computing speed. The signal 
conversion process is shown in figure 5.7. 


The references of pathology images.were.extracted from real signals which are 
shown in figure 5.6 proves how easy it can be identified with the original 
pathology. Of course, the graphics are enhanced by the computer through the 
interpolation procedure where pairs of neighboring points are plotted using anti- 
aliased filter to reduce the pixelization. In fact, without being processed the 
signals are the same as in the matrix shown in figure 5.7. The construction of the 
above shown algorithm is achieved through a very simple calculation routine that 
consists only of a few lines of program. Reference signals were rescaled to | Volt 
for a total of 50 points. The x axis number of points remains unchanged against 
the original signal and dynamically allocated between 100 and maximum 350 
points. Under these conditions the reference image is represented in a matrix of 
50x350 points. To test the recognition algorithm, the routine uses computer 
generated pixels of the image without intensity parameter Z (pixel has value 1) 
and no interpolation therefore a minimum number of points. The lowest number 
of points (42) is contained by the reference of normal heart rhythm. Analogous to 
the procedure for the reference signal, it applies the same above procedure for the 
analyzed signal. As a pattern recognition method, we have used the computation 
of the Euclidian distance between the analyzed signal and the stored references. 
The Euclidian distance is computed using the equation 5.5: 


d(p,.4;)= 3 >, -p,} 


ml al (5.5) 


6. Conclusions 


It was presented a set of algorithms for processing of the PCG, which are part of 
two main methods of pathologic classification, along with the obtained results. 
The first method involved the analysis of spectral content of the PCG signal and 
then based on a complex algorithm; the frequency spectrum is displayed on a 
color LCD display. This method of evaluation of the signal is addressed to 
medical personnel who can understand the time-frequency scalogram created. The 
PCG signal correlation method includes an algorithm for accurate heart rate 
detection using autocorrelation function. Determining the precise heart rate is 
dictated by the need of temporal rescaling of the wavelet decomposition at level 
three, to bring the acquired signal to the same HBR value of the reference signal. 
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The last step is to precisely detect the envelope of the PCG signal by using the 
discrete signal Hilbert transform. The PCG signal correlation is prepared by 
converting the signal into a two-dimensional vector, making the actual signal in a 
two-dimensional image. This method opened the way to the methods of imaging 
processing especially for shape recognition. The sum of these presented numerical 
processing methods is an important step in the processing of a signal virtually 
unapproachable by traditional methods of numerical processing. 
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